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FOREWORD 


This report was prepared by UNIDEV , INC. for the 
George C. Marshall Space Flight Center, National Aero- 
nautics and Space Administration, under Contract NAS8-21418 
entitled "Wind Induced Oscillation of the Saturn V Vehicle 
Caused by Realistic Ground Winds." The work was conducted 
under the technical direction of the Aero-Astrodynamics 
Laboratory, with Mr. G. A. Wilhold acting as technical 
monitor . 

The study was executed by the Solid and Fluid 
Mechanics Department, directed by Dr. L. L. Fontenot. 

Mr. T. Coffin served as Principal Investigator. Computer 
Program development was performed by Mr. T. S. Chandler. 

This document deals with the computational 
aspects of the "Ground Winds" problem. Theoretical con- 
siderations are presented in reference /!/. 
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ABSTRACT 


This document describes a Digital Computer Pro- 
gram developed to estimate Saturn V response to pre-launch 
wind conditions at Cape Kennedy. The program computes 
displacement and bending moment statistics as a function 
of parameters defining the atmospheric environment. A 
sample problem is provided to illustrate utilization of 
the program. 
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INTRODUCTION 


An analysis of space vehicle response to pre- 
launch winds at Cape Kennedy is presented in /l/. The 
vehicle is assumed to be axisymmetric and supported 
symmetrically. The equations of motion are derived 
based on simple "strip theory." The aerodynamic forcing 
function is assumed representative of a non-homogeneous , 
stationary normal process. Generalized harmonic analysis 
is applied to obtain response statistics. Lift and drag 
response are computed separately. Joint statistics are 
then computed to obtain an estimate of the probability 
density and cumulative distribution of resultant (i.e., 
vector sum) displacement and bending moment. 

This document describes utilization of the com- 
puter program resulting from the analysis. A sample 
problem is included to illustrate application of the 
program. 
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SECTION 1 


SUMMARY OF SIMULATED EQUATIONS 
An analysis of Saturn V response to pre-launch 
wind conditions is documented in reference /l/. The 
equations programmed for computer simulation are sum- 
marized in this section. 

1.1 NOTATION 

x Coordinate indicating position along vehicle 

length (identical with station no.). 

y Displacement coordinate in the drag direction 

(parallel to mean wind vector). 

g Displacement coordinate in the lift direction. 

\p n Displacement mode vector. 

y n Bending moment mode vector. 

M Generalized mass, 

n 

8 Generalized structural damping ratio, assuming 

n viscous modal damping. 

v n Undamped natural frequency. 

D Mean diameter of discretized vehicle model, 

p Density of atmosphere. 

C^ Drag coefficient. 

R Reynolds no., based on local mean wind speed and 

vehicle diameter. 

V Mean wind speed. 

A Peak reference wind speed at the 18.3 meter 

level from grade . 

k Empirical statistic defining the wind speed 

profile at Cape Kennedy. 
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V 

Frequency . 

$ yy ’ ^zz 

Cross spectral density of wind turbulence in 
drag and lift directions, respectively. 

D ,M 

y y 

Mean drag displacement and bending moment, 
respectively . 

‘f’yy ’ ‘f’ym 

Mean-square density spectrum of drag displace 
ment and bending moment, respectively. 

2 2 
a , a 

y ym 

Variance of drag displacement and bending 
moment, respectively. 

2 2 
a , a 
z ’ zm 

Variance of lift displacement and bending 
moment, respectively. 

l|D||,||M|| 

Maximum absolute displacement and bending 
moment, respectively, associated with joint 
occurrence of 3a response. 

r 

Absolute (resultant) displacement or bending 
moment, as indicated. 

f (r ) 

Probability density of absolute displacement 
or bending moment, as indicated. 

FCr) 

Cumulative distribution of absolute displace- 
ment or bending moment, as indicated. 


A number of the above quantities are not specified, but 
implied, since highly redundant equations for displacement 
and bending moment statistics are not repeated. Coding 
terminology for required input is defined in the section 
on computer program utilization. 

1.2 DRAG RESPONSE 

Mean drag displacement is computed by 
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where 


— _ on oo« /-0.025X + 28.96 ^k < n A 
V(x) = 39 . 37A,. { ( 18^ ' \ 1 + 


f _ 18.3 v 

'•0.025X + 28 .96 ; 


(0.283 - 0.435 e 0,2A * ) | 1 


1.98 - 1.887 e 


-0.2A, 


(1-A) 


and 


1.0 


,R(x)<2 . 8 x 10' 


c d (x) 


3.155(5. 7 624-log 1Q R(x) ) , 2 . 8 x 10 <R(x)<5 x 10' 
0.514(log 10 R(x)-5.3099) ,5 x 10 5 <R(x)<3 x 10 6 


0.6 


,R(x) >3 x 10 


(1-B) 


Mean drag bending moment is defined by a similar expression’ 


py (x) 

M (x) = > -~-2 - 

y Z-J 8tt 2 v 2 M 0 

n n n £ 


f 


C d (x)D(x)V 2 (xH n (x)dx 


( 2 ) 


The variance of drag displacement is estimated by 


a 2 (x) 

y 


J <f> yy (x;v)dv 


(3) 


where 


<j> (x . ; v ) 

yy i 
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p2 ^n^ X i^ I | [4tt 2 v 2 (M^+ / D 4- (x)^(x)dx) ]' 


/' 
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,2 , . ,2 

/ u. • - i l > (x)ii 1 

n n 4 J r n 
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[ (l-(— ) 2 ) 2 + (■ 


M 


n 


n 


M 


n + D 2 (x)i|/ 2 (x)dx 

£ 


> (6 n + 


n 


3 



V / c , (x)D(x)V(x)^ 2 (x)dx) 2 (— ) 2 ] 


TT V M 

n n 


// 


X II I C d^ X l^ C d^ X 2^ D ^ X l^ D ^ x 2^ V ^ X l^ V ^ X 2 ^ 


+ (^p- D(x^)D(x2 ^ ^ 


$ yy (x l 5 x 2 ;v)^ n (x 1 )^ n (x 2 )dx 1 dx 2 


The cross spectral density of wind turbulence is defined in 
references / 2 , 3 / . In the notation of these references: 


OyyCx-j^ ,x 2 ;v) = 1.55 x 10 S(h^,h 2 ;v) 


h. = 0.025x. + 28.96 

l l 


Now the (real part of) cross spectral density is: 


-0.3465 


S(h 1 ,h 2 ;v) = ^S(h 1 ;v)S(h 2 ; v) e 


cos ( 2TryAf ) 


vh 2 vh^ 

W(h 2 ) WO^) 


Af n n = 0.036 

U • o 
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Y 


0.7 


0.3 


h l +h 2 


h l +h 2 


-$100 
->10 0 


and 


W(h) = A* ( 


<P^> 


(0.283 - 0.435 e °* 2A *> "j _1 


*'18.3 


) J 


1 + 


1.98 - 1.887 e - 0 ‘ 2 A * 


( .3-D ) 


The local spectral density function required is represented 
by 


S(h;v) 


BaU^C, f/f 
“ - 1 m 

v(l + 1.5C 0 (f/f ) C 2 ) fe/3C 2 
2 m 


( 3-E ) 


where 


f = 


vh 

W(h) 


h 

f = C a ( T g- T ) 
m 318.3 


6 - ( Jl — )° 5 

B * ( 18 . 3 ; 


u* = C g V?(18.3) 
C 1 = 6.198 
C 2 = 0.845 
C 3 = 0.03 

c 4 = 1.00 

C 5 = - 0.63 

Cg = .0886 
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The variance of drag bending moment is similarly 

computed : 

00 

o 2 = /<f>(v) d v (4) 

ym / Y ym 

J 0 

Computation of <)> is identified to that for $ , with 

2 2 
ip (x.) outside the brackets in (3-A) replaced by y (x.). 
n l n l 

1.3 LIFT RESPONSE 

The variance of lift displacement is calculated 
in the same manner as (3) with 


<f> (x . 
y zz 


^;v) = p 2 ^ ( x.) | j [4rr 2 v 2 (M n + jp- J” D 2 (x)i(j 2 (x)dx) ] 2 

n^ ( ' Z 


C ( l -<— ) 2 ) 2 + ( 
v n M 


M 


n 


n+ it y D (x)*. n < x 


-) (e„ + 


2 , 2 /-jdx '^ n 


/ 


-1 


4tt v M / C d (x)D(x)V(x)K; 2 (x)dx) 2 <£-> 3 

n n / n 


x J ff | c d < x 1 > c d ( x 2 )D(x l )D(x 2 )V(x i )V(x 2 )+ J (v7rD(x i )D(x 2 )): 


$ zz (x l ,x 2 ’ v) ^n (x l^n (x 2 )dx i dx 2^ 


(5) 


and the following constants defining the cross spectral 
density of turbulence in the lift direction: 


Af q g = 0.045 


6 



1.4 


V h 2 


2 * 


<100 


Y = 


h + h 

0.5 , " o " >100 


C 1 = 3.954 
C 2 = 0.512 

c 3 = 0.1 

C 4 = 0.58 
C 5 = 0.35 
C_ = 0.0886 

b 


The variance of lift bending moment is computed similar to 

2 

displacement, with y n ) inserted in (5) in place of 

2 

\jj n (x^), outside the brackets only. 


1.4 RESULTANT RESPONSE 


Maximum absolute displacement associated with con- 
current 3a lift and drag response at a given location is 
computed by 


D 


D 2 


9 + 


2 2 
a -a 
z y 


1/2 


I D | a 

» — 9- <1> and a >0 

3 <°H> y 


C | D | + 3 CTy ) » otherwise 


( 6 ) 


Maximum absolute bending moment associated with 

concurrent 3a lift and drag response at a given location 

is obtained by (6) with the substitution 

D -►FT 

y y 

a —►a 

y ym 

a —►a 

z zm 
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The probability density of absolute (resultant) 


displacement at a given vehicle location is defined by 
2 7T 


f(r) 


2 it a a 

y z o 


/ 


e Q d0 


r> 0 


(7) 


= 0 


0 


where 



Integration of (7) yields the cumulative distribution of 
absolute displacement: 


F(r) 



f ( r ' ) dr ’ 


( 8 ) 


Equations (7) and (8) are used to define the same statistics 
associated with absolute bending moment with the substitution 



a —►a 

y ym 


a —►a 

z zm 

It should be noted that (7) and (8) represent conditional 
density and distribution functions, respectively, based on 
the assumed values taken on by meteorological constants. 

In the following sections, organization and 
implementation of the computer program developed to solve 
the above equations are described. 
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SECTION 2 


PROGRAM DESCRIPTION 

Derivation of the equations "leading to the com- 
puter program described herein is documented in reference 
/l/. The wind turbulence experienced by the vehicle is 
assumed representative of an inhomogeneous, stationary 
normal process. In the computation of joint response 
statistics , lift and drag components are assumed inde- 
pendent . 

The numerical integrations are carried out by 
application of Simpson's method and in some cases by the 
trapezoidal rule. The program offers the user several 
options. By setting a flag card, analysis of drag response, 
lift response, or both, including combined response, may 
be obtained. Stacking of runs is provided for, permitting 
a series of analyses to be performed with variation of 
any combination of input parameters. 

Printout includes a definition and listing of 
input terms and computed results. Results are also written 
on tape and processed by the SC-4020 plotter. The program 
was written and checked out on the UNIVAC 1108, Exec. 8 
system. 
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SECTION 3 


PROGRAM UTILIZATION 

3.1 GENERAL PROGRAM ORGANIZATION 

Organization of the computer program, input data 
and case stacking scheme is illustrated schematically as 
follows : 

| Control Cards 

{ Input characters to be used as an identifying 
title Cuse columns 7-72) 

| Main Program and Subroutine 

CASE DATA 

Second 

Column 

$NAM1 

IFLAG=3 

(Complete set of input data - first case) 

$ 

$NAM1 

(Load only parameters to be changed for second case) 

$ 


IFLAG = 0 } St °P s 

$ 

| Control Cards 

In the above iterative scheme any number of parameters may 
be varied for each case. As explained in Section 3.2.1 
"IFLAG=3" requests an analysis in lift and drag directions 
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and the computation of resultant (joint) response statistics. 
A description of the coding format for a sample problem is 
provided in Section 3.2.4. 

3.2 DESCRIPTION OF INPUT 

Three sets of input data are required which define 
general program instructions, meteorological parameters and 
structural dynamic characteristics. The program is designed 
to read input in Namelist format. Therefore, the exact allo- 
cation of all spaces on a coding form is not required. 
Definition of constants must be separated by a comma and 
the elements of arrays must be separated by commas. Data 
should be entered beginning with column two of each line on 
the coding form, i.e., column one must be left blank. 

3.2.1 General Instructions 

The instructions define general requirements of 
the problem to be solved. All units are in the pound-inch- 
second system. 

! = 0 Stops program after last case 

= 1 Compute only drag response 

= 2 Compute only lift response 

= 3 Compute drag, lift and resultant response 

NIPD Number of vehicle locations for which resultant 

displacement statistics are to be computed (zero 
unless IFLAG=3). 

IPTD Station number (in hundreds) of locations at 

which resultant displacement statistics are to 
be computed (zero unless IFLAG=3). 
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NIPM 


IPTM 


SBARD 


SIGDS 


SBARM 


SIGMS 


Number of vehicle locations for which resultant 
moment statistics are to be computed (zero 
unless IFLAG=3). 

Station number (in hundreds) of locations at 
which resultant moment statistics are to be 
computed (zero unless IFLAG=3). 

Mean lift displacement due to vortex shedding 
from scaled wind tunnel data at points IPTD 
(zero unless IFLAG=3). 

Standard deviation of lift displacement due to 
vortex shedding from scaled wind tunnel data at 
points IPTD (zero unless IFLAG=3). 

— 6 

Mean lift moment (lb.in.X10 ) due to vortex 
shedding from scaled wind tunnel data at points 
IPTM (zero unless IFLAG=3). 

Standard deviation of lift bending moment 
(lb . in .X 10 - 6 ) d U e to vortex shedding from 
scaled wind tunnel data at points IPTM 
(zero unless IFLAG=3). 


3.2.2 Definition of Meteorological Conditions 

The quasi-static and turbulent wind distribution 

at Cape Kennedy is characterized by two input constants. 

XLAM Peak wind speed at the 18.3 meter level 

(meters/second) . 

XK A constant defining the wind speed gradient. 

Values of the above statistics for varying levels 
of risk and exposure time are available from the M.S.F.C. 
Aerospace Environment Division. See reference /4/. 


3.2.3 Structural Dynamic Characteristics 

The launch vehicle to be analyzed in a given 
propellant loading condition is synthesized by the fol- 
lowing parameters. All units are in the pound-inch-second 
system. 
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N Number of structural' modes to be considered. 

D Mean diameter of discretized vehicle segments. 

This array must be entered in the order of 
increasing vehicle station number. 

XNUN Undamped natural frequencies (Hz). These must 

be entered in the order of increasing mode 
number . 

2 

XMN Generalized masses (lb. sec. /in.). These must 

be entered in the order of increasing mode number. 

BETAN Generalized structural damping ratios. These 

must be entered in the order of increasing mode 
number . 

PSI Displacement mode vectors (in./in.), normalized 

to tip displacement. 

XM Bending moment mode vectors ( ^£ — *X 10~^ ) , 

normalized to tip displacement. 

The constants PSI and XM obviously constitute two 
dimensional arrays. PSI(i,j) represents the displacement of 
element i in mode j. The same convention applies to XM. 

For a given mode number the modal values must be entered 
in the order of increasing segment number (station number). 

3.2.4 Input Format 

The coding format for a sample problem is illus- 
trated by the following printout. The input represents a 
series of four analyses with a fixed reference wind speed 
but varying wind speed gradient. Note that no information 
is entered in column 1. In repition, the input data is in 
Namelist format so the exact allocation of fields for the 
entry of constants is not necessary. For example, consider 
the fourth row of input. These four constants may be 
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entered in any order and in any column other than the first . 
These same four constants may also be entered in four 
separate rows (instead of one), yielding a separate data 
card for each constant. The Namelist data format was 
selected to provide flexibility to the program user. 


14 



SNAM1 
I FLAG=3 
N=3 

I PTP=i ,IPTD*42,MPP*l,MFD=l 
XK.,278 ,XLAP=10.3 
XMN*3. 273 ,0.782, 1.392 
XPN =780 .4,127.52, 37 .41 
BETAN=3*0.01 

0=713. ,594. ,436. ,13*428., 396., 8*396. ,347. 8, 286.8 , 297 ., 297. , 2 59. 8 , 25 9. 8 , 
259.8,247.5,217.8,182.2, 154., 154., 128 .7, 3*26 ., 5 . 

PSI( 1,11=9.9 ,0.01 ,0.02, 0.03, 0.05, 0.05, 0.06, 0.08, 0.1, 0.12, 0.14,C. 15, 

0.17,0.2 ,3.21 , 0 . 24, 0 . 25,0. 27, 0. 30, 0.33, 0.35,0. 37, 0.40, 0.42, 0.45,0.4 7, 

0. 50,0. 53, 0.55,0. 58,0. 60,9.64,0. 66, 0.7, 0.73, 0.75, 
0.80,0.84,0.86,0.90,0.95,0.99 

PSI ( 1, 2) =-0.0 1,-0 .02, -0.33, -0.045, -0.05, -0.055,-0. 66, -0. 67, 6*-0.08 ,-C .U75 , 
-0.07,-0.065,-0. §6,-0 .05,-0 .04, -0.03, -0 .015,-0. CCS, §.C1»C. 025, C. €5 , 
0.07, 0.09, 0.12,0. 15,0.155,0. 2, 0.25, 0.2 59, 0.345, €.4, £.4 7 5,S.55, £.65, 
0.748 ,0.85 ,9 .96 

PSI (1,31=0.01,0.02 ,0.025, 0.03, 0.04, 0.045, 3*0 .048,0 «C 4,0.03 »C .025 ,0 .£2 , 
0.01,0.0,-0.015,— 0 .02,-0 .03,-0 .04,-0 .05,-0.055,-0.06 ,-C. 075 ,-C . 09 , 
-0.095 ,5* -0. 1, -0 . 0 8, -0. 0 7, -0 .05, -0.C2, 0.0 2,0 .09,0.15,0.25, 0.38, C. 55 , 

0 .75,0.95 

XM1 ,11* -13. 4 ,-12. 5, -11. 7, -11. 0,-10. 1,-9 .4, -8. 7, -8. 0,-7. 2,-6. 5, -5. 8, -5. 3 , 
-4. 7, -4. 2,-3 . 75,-3. 4, -3. 0, -2. 6,-2. 25, -2. 05,-1. 8, -1. 65,-1. 45, -1.25, 
-1.10 ,-0 .85 , -0.70,-0.60, -0.50, -0.3 5,-0 .27,-0. 24,-0. 2,-0. 15,-C.l » 
-0.07,6*0.0 

X P ( 1 » 2 1 = 7.84 ,6.0,4.45,2.60, 1.0, -0.05,- 1.9, -3. 20, -4. 55, -5. 7, -6. 4,-7. C, 

-7. 5, -7. 65, -7. 5, -7. 35, -7. 24, -7. 1,-6. 75, -6. 3,-5. 9, -5. 4, -5. 0,-4. 5 ,-4.06, 
-3. 7, -3. 2, -2. 75, -2. 4, -2. 05, -1.75, -1.45, -1.18,-0.90, -0.7, -0.45, -G. 23, 
-9. 14 ,-0.10,-0.02 ,2 <0.00 

XP( 1,3) *-3.65, -2. 0,-0.24,1.20, 2. 63, 3.80, 4. 59, 5. 40, 6. 20, 6. 42 ,6.20, 5. 77, 

5. 22 ,4. 39, 3.40, 2. 45, 1.45, 0.60, -.18 ,-0.6, -1.0 6,-1. 40,-1. 80,-2. 1C, 

-2. 39, -2. 60, -2. 80, -2. 95, -2. 80, -2^60, -2. 35, -2. 08, -1.83, -1.45, -1.2% 

-0. 82, - 0 . 50, -0. 30, -C5. 20, -0.15, 0.04, 0.00 
SIGNS (1 1=0.0, SEARP(1)*O.0,SIGCS 111 * 0 . 0,SeARC( 11 * 0.0 
t 

SNAN1 
XK=. 216 
t 

SNAM1 
XK=, 153 
$ 

SNAM1 

XK=.09O 

$ 

SNA Ml 
I FLAG*0 
S 
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3.2.5 Additional Constants 

A number of constants are utilized which generally 
remain fixed for a parameter study concerning, say, various 
wind speeds or propellant loadings. In order to minimize 
required input, these constants are included in the main 
program. 

NODE Number of elements in discretized vehicle model. 

XI Station number of first element from fixed end 

of vehicle Cin.). 

DX Distance (fixed) between elements (in.). 

2-4 

RHO Mass density of atmosphere (lb. sec. in. ). 

The presently programmed constants represent a 
42 element vehicle model with elements centered at stations 
100 ,200 ,... ,4200 . This model may be altered by modifying 
the above constants in the main program. Consistency with 
input data and program dimensions must, of course, be 
maintained . 

3.2.6 Program Restrictions 

The program is dimensioned for a maximum of: 

4 bending modes, 

42 vehicle elements, 

10 vehicle stations at which resultant proba- 
bility density and cumulative distribution 
functions are to be computed (for displace- 
ment and bending moment, separately). 

The lift and drag spectra of turbulence are 

computed based on /2,3/, with constants appropriate to 

neutral wind conditions. 
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3.3 


DESCRIPTION OF OUTPUT 


The program output consists of the following 
printed information: 

A. Case Information 

1. Program title. 

2. Definition of input terms and appropriate 

units . 

3. Listing of input values. 

B. Results 

1. Definition of output terms and units. 

2. Output in drag direction: 

a) Mean displacement. 

b) Mean bending moment. 

c) Displacement standard deviation and 
variance . 

d) Bending moment standard deviation 
and variance. 

e) Mean wind speed distribution. 

3. Output in lift direction: 

a) Displacement standard deviation and 
variance . 

b) Bending moment standard deviation 
and variance. 

4. Resultant response 

a) Maximum absolute displacement asso- 
ciated with 3 sigma lift and drag 
response . 

b) Maximum absolute bending moment 
associated with 3 sigma lift and 
drag response. 

c) Probability density and cumulative 
distribution of displacement at 
vehicle locations specified. 

d) Probability density and cumulative 
distribution of absolute bending 
moment at vehicle locations specified. 
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The functions (c) and (d) above will generally 
be computed for differing vehicle locations. When joint 
statistics are requested at a number of points, all dis- 
placement statistics are printed out in the order of dis- 
placement locations specified, followed by bending moment 
statistics in the order specified. 

Computed results are also written on tape and 
processed by the SC-4020 plotter. The order of data 
presentation is the same as described above. 
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SECTION 4 


SAMPLE PROBLEM DESCRIPTION 

A sample problem is presented to illustrate com- 
puter program operation. Consider the format defined in 
Section 3.2.4. This format indicates: 

1. Structural dynamic characteristics representing 
the Saturn V (SA-503) vehicle in the fully 
fueled condition, including three bending modes. 

2. A reference peak wind speed of 10.3 m/sec. 

Cat the 18.3 m. level), which is equivalent 
to a mean wind speed of 20 knots. 

3. Lift, drag and resultant response are to be 
computed . 

4. Resultant displacement statistics are to be 
computed at station 4200. 

5. Resultant bending moment statistics are to be 
computed at station 100. 

6. Four cases are called for, with the value of XK 
being varied and all other parameters held 
fixed . 

As stated previously, the quantity XK is an 
empirically determined statistic which defines the wind 
speed profile shape at Cape Kennedy, see reference /4/ . 

The values selected for this example represent vehicle 
response to the 3a, 2a, la and mean wind profile envelope, 
respectively, for a fixed reference wind speed. 
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SECTION 5 


SAMPLE INPUT 

The printout of input data is entirely redundant 
for all cases, with the exception of XK, which is specified 
in Section 3.2.4. Therefore, only the first case is pre- 


sented . 



number of points at which joint probability is desined cdisplacementi 
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SECTION 7 
PLOTTED RESULTS 

A series of 12 response quantities are plotted 
for each case, yielding 48 figures. Since results for the 
four cases are printed in Section 6, only the first case 
is illustrated. 
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SECTION 8 
COMPUTER LISTING 
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SATURN 5 RESPONSE T2 FRE-LAUNCH WINC L0AOS 

VEHICLE SA-503 

FULLY FIELEC C0NCITJ2K 


DIPENSI2N ABSC(42),AESM(42) 

DI PENS I 2 N H(42)»W(42)tV(42) ,V2 ( 4 2 ) , FM( 42 ),B£TA( 42) ,0 (42 I , 

IX(42) ,CD (42) ,C2(42) ,CDV(42 )» CCV2<42)»XNLN< 4),XMN< 4), 

«BETAN(4) ,PSI (42.4 ) , FS 12 ( 4 f 42 ) , C0N5( 4 ,4 2 ) ,C0N 6 ( 4 , 42 ) , 

$ T A B 1 ( 4 2 ) * T A 8 2 ( 4 2 ) t S l M 1 ( 4),SUM2( 4),C0N1( 4), 

*C0N3( 4) ,XNUa03),S(42),PH(4 2»42)tPHIYC 1€0), 

$ARPAY(42) ,CB(42) , S IC CZ ( 4 2 ) * S IGDY ( 42 ) , $ IGMZ ( 4 2 ) , SIGMY (42 ) ,. 

*XNU2 (100 ) »SUP3( 100*4 ) , FF (4 2 ) , CCINTf 2 ) , D INTC 150) 

DI PENS I 2 N C0VS(42)*C4(42 ) » SUM4 ( 10C»4 ) » h 182 ( 42 ) 

DI PENS I 2 N C2N7(4»42)»XP(42*4 )»XMB(42)» 

$C0N8(4,42> ,PHIM( ISO ),XP2 (4*42 ),SFANS( 150, FANS ( 150, 

*SIDE9(12 ) . B0TTP ( 12 » 2 ) , IFTM( 1C )* IPTD( 10 ) 

DIPENSI2N SICEUU2) 

DIPENSI2N AC ARY (6) ,VARC(42 ) ,V ARM ( 42 ) , B0TT( 12 ) t S IDE1 ( 12 ) • SIDE 2 ( 12 ) ♦ 
$SICE3 (12) ,SICE4(12) tS ICE5( 12 ),SICE6( 12 ) f S IDE 1 ( 12 ) * SIDE 8 ( 12 ) 

DI PENS I 2 N e2TP(12)»E0TPP(12) 

DI PENS I 2N SIGPS142)»SEARP(42)»SIGCS(42)»SBARC(42) 

ECLIVALENCE ( e0TT P ( 1 , 2 ) , E0TPP ) 

ECUIVALENCE (B0TTP ( I* 1 ) , E0TP ) 

DATA ACARY/36H 401240 E IN 222 T CHANDLER APR 70/ 

DATA S I C E1/3*6H ,36FPEAN ORAG DISPLACEMENT - INCHES * 

*3 *6H / 

DATA SICE2/2«6H ,42FMEAN DRAG BENDING MOMENT LB« IN. X 1C-6 

* ,3*6H / 

DATA SI C E3/3 *6H ,3fcFVARIANCE 0F CRAG DISPLACEMENT IN.2 .* 

$3 *6H / 

DATA SIDF4/3*6H ,36FVAR|ANCE 0F LIFT DISPLACEMENT IN.2 * . . 

$3*6H / 

DATA SICE5/6H ,54HVARIANCE 0F CRAG BENDING M0MENT LB. 2 IN.2 X 

$ 10-12 * 2 «6H / 

DATA SICE6/6H ,54FVARIANCE 0F LIFT BENDING M0MENT LB. 2 IN.2 X 

$ 10-12 ♦ 2 «6H / 

DATA SI CE7/3 »6H ,34MeS2LUTE DISPLACEMENT INCHES ,4*6H 

$/ 

DATA SICE8/3*6H ,42HAES0LUTE BENDING M0MENT LB. IN. X 1C-6.. 

$ »2*6H / 

DATA SICE9/4«6H ,30F CUMULATIVE C ISTR I8UTI0N *3*6H 

%/ 

DATA SICE10/5«6H » 24HPR2EABIL LTY DENSITY *3*€H / 

DATA B2TT/4«6H ,24FVEHICLE STATI0N - INCHES, 4*6H / 

DATA BeTF/4*6H ,3OHABS0LUTE DISPLACEMENT INCHES ,3*6H .... ... i 

DATA B2TFP/2-6H »42HABS£LUTE BENDING M0MENT LB. IN. X 1C-6 

% ,3*6H / 

STANT=1 . 

WRITE (6,6051 ) 

6051 F0RPAT (IHl »4X42HSATURN 5 RESPeNSE T0 PRE-LAUNCH WIND L0AOS//) 

WRITE (6, 6061) ... 

6061 F 0PPAT ( 4X19HCEF IN IT I£N 2F INPUT//) 

WRITE (6, 6001) 

6001 F0PMAT (4X25HN = NUMBER EF M0OES) 

WRITE (6, 6006) . . 

6006 F0RPAT (4X53HMPM = NUMBER 0F P0INTS AT WHICH J0INT PR0BABILITY, 

II X19HIS DESIRED ( M6MENT ) ) .. ... . 

WRITE (6 »6025 ) 

6025 F0RMAT (4X53HNIPC = NUMBER £F P0 INTS AT WHICH J0INT PR0BABILIIY* 
11X25HIS CESIRED ( C ISPL AC EM ENT ) ) 6Q 



WRITE (6,6008 ) 

6?C8 F2RMAT (4X46HIFTM = FEINTS AT WUCH J2INT PR0B ABILITY IS, 

II X16HCESI REC (MEMEM ) ) 

WRITE (6 ,6C09 ) 

6 ■ j 9 FERMAT (4X46H IPTC = POINTS AT WHICH J0INT PROBABILITY IS, 

11 X22HCE S I RE C ( C I S P L ACEH ENT )) 

WRITE (6 ,6010 ) 

6'l .i FERMAT (4X37HC = AVERAGE CIAMETER 0F VEHICLE, 

11X13HSFGMENTS-- IN. ) 

WR ! TE (6 ,6(333 ) 

603 3 FERMAT ( 4 X 4 7 H X K = CCNSTANT DEFINING WIND SPEED GRADIENT) 

WRITE (6,6020 ) 

6520 FERMAT (4X52HXLAM = PEAK WINC SPEED AT 18.3 METER LEVEL — M/SEC) 
WRITE(6,6021) 

6*21 FERMAT (4X42HXNUN = UNCAMFEC NATURAL FREOLENC I E S — HZ) 

WRITE (6 ,6022 ) 

6522 FERMAT (4X43HXMN = GENERALIZED MASSES— LB*SEC**2/IN.) 

WRITE (6 ,6023 ) 

6523 FERMAT (4X35HEET AN = S'! PUCTLIR AL CAMPING RATIOS) 

WRITE (6,6024) 

6124 F 0 P M AT (4X43HFSI = DISPLACEMENT M0CE VECTOR — IN. /IN.) 

WRITE (6 ,603 0 ) 

6130 FePNAT (4X50HXM = FENCING MOMENT M0CE VECT0R — ( LB* IN . / I N. ) , 

14H*E-6) 

WR I TE ( 6 , 60 26 ) 

602 6 FORMAT (4X43FSIGMS =? STANCAPC CEVIATI0N 0F LIFT MOMENT , 

11 X34HCUE Te VORTEX SH E C C INC — L B* IN .*E- 6 ) 

WRITE(6,6027) 

6027 FEPMAT (4 X49HSEARM = MEAN LIFT MOMENT CUE T0 V0RTEX SHEDOING, 
llX12H--Le*IN.«E-6) 

WRITE (6 ,6028 ) 

6028 FORMAT (4X49HSIGCS = STANCAPC DEV I ATI0N 0F LIFT DISPLACEMENT, 

1 1 X27HDUE Te VERTEX SH E C C ING-- IN . ) 

WRITE (6,6029) 

6 029 F0RMAT (4X39HSBARC = MEAN LIFT DISPLACEMENT DUE T0 , 

11 X20HVePTEX SHEDDING— IN.) 

N0CE = 42 
NFREC =85 
XR = 3. 

NFD = 100 . . 

DTHETA= .62851853E-1 
RHe = • 1 1470000 E -6 

XI = 150. 

DX = 100. 

FREC1 = .052 

DFREC = .017 

3500 F0PMAT ( 1 H 1 » 1 7 H N NIPM MPD/3I5, 

1 / /4X ,6H I FT M =,(8115)) 

3001 FERMAT (/8X2HXK, 14X,4HXLAM/2E16.8) 

3CC2 FERMAT (/23X,4HXNUN/5E16.8) 

3003 FERMAT (/23X ,3HXMN/5E 16 . 8 ) 

330* FERMAT (/23X , 5HBET AN/5 E 16 .8 ) 

3005 FERMAT (/50X,1HD/(6E16.8 ) ) 

3*306 FERMAT (/5QX»3FiPSI/(fcE16.8) ) 

3007 FERMAT (/50X,2HXM/ (6E 16.8 ) ) 

3008 FERMAT (//4X,6HIFTC =,(8115)) 

3009 FERMAT(3X,7HSIGMS =,(8E15.7)) 

3010 FERMAT (3X,7HSEARM =,(8E15.7)) 

3011 FEPMAT(3>,7HSIGCS =,(8E15.7)) 

3012 FERMAT (3X,7HSEARC *,I8E15.73) 

CALL I CE NT ( ACARY *9 ) 
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NAPELIST/NAVl/N, NIIPM,NIPD, IPTM,IPTD,D, 

tIFLAG,XK , XLAM,XNLN,XMN,RETAN,PSI,XM 

«, SIGN'S, SPARP,SIGCS,S8ARC 
4 RE AC ( 5 » N A M ) 

I F L AG = 1 V CIPECTIPN, IFLAC=2 Z CIRECTI0N, IFLAG=3 B0TH DIRECTIONS 
IFLAG=:, LAST CASE 
I B PTH = 

1 IF (IBPTF )3 ,3 ,50.' 

3 IF (IFLAC )813,815,2 

2 CONTINUE 

WRITE (6,3^;.!-i) N y NIPM,NIPD, 

$ (IFTP(I) ,I = 1,MFP) 

WRITE (6,3009 MS IGPS ( I ) , 1 = 1 , MPM ) 

WRITE(6»301 MSEARNMI ) » I = I,NIFM) 

WRITE (6,3008 ) ( I PTC ( I ) , 1 = 1 ,NIPC ) 

WRITE(6,3011)(SICCS(I),I = 1,MPC) 

WRITE(6,3012MSBAPC(I),I = 1,MPC) 

WRITE (6 ,30.1) XK, XL AM 

WRITE (6 ,3032) ( X NU N ( I ) , I = 1 , N ) 

WRITE (6 ,30)3) (XPN(I ) , 1 = 1, N ) 

WRITE (6 ,3 004) ( E ET AN ( I ) , 1= 1 , N ) 

WRITE (6 ,30:5) (C(I),I=1,N2CE) 

DP 6 1=1 , N 

WRITE (6 ,3006 ) ( FS I ( J » I > » J= 1 » N 2 C£ ) 

6 WRITE (6 ,3007 ) ( X M J , I ) , J= 1 , NP C E ) 

X2=X1+N2CE«CX 
FREC11=FPEC1-CFREC 
De 10 1=1 » NF PE G 
X I = I 

XNL (I )=FPECi 1+XI*CFF£C 
10 XNL2(I)=XNU( I ) * « 2 
N1=N2CE-1 
D X 2 = D X «C X 
PI =3. 14159265 
P l 22 -P l * *2 

P1=4.0«FI22 
P2 =2 • ; * F I 
P3 = FI/4. 

P4 = RHE«F 1/4. .0 

NPT= (PI/C THETA ) 4 1.00301 

RH22=PH2«RH2 

P6 = PH£/ ( 8 • 0 * F 1 2 2 ) 

ELAP=EXF ( -0 . 2 « X L A P ) 

P = C.283- ;.435«El AM 
GJ=1.98-1.887«ELAP 
W1P3=XLAP/(1 ,0+1 . 0 / GO ) 

X11=X1-CX 
De 50 I=1,NPCE 
X I = I 

X ( I )=X11+XI«CX 
H(I)=0. 25«x { I ) +28.956 
HI 83 ( I ) =F ( I ) /18 .3 

W(I)=XLAN*H183( I )«*XK/( 1.0 +( 18,3/HI I) )«*P/GC ) 

V(I) = 39.37«Wm 
V2 (I )=V ( I ) *V ( I ) 

R = 44.2321*V( I)*C( I ) 

IF (P-2.8E5) 20,20,25 
20 CD (I)* 1.0 
Ge T £ 46 

25 IF ( P-5 .0 E 5 ) 30,30,35 

3C CD (I)=3.155»(5.76 24-ALeG10(P ) ) 62 



r>£ T? 4 6 

i5 IF;F-?. 6 ) 4, ,45,4>, 

4j CD <J )= '.514MALgGlC(R)-5.309S) 
ce re 46 
a c c r ( i ) = : . 6 

46 CeMIfsLt 

0 2 ( I ) = C ( I ) * C ( I ) 

04 (I )=D?(I ) «C2 ( I > 

rcvuj*ccm«c( i)*vu ? 

CDVS(I)=CCV(I)*«2 

c o v 2 ( i )».ccvm«v<i ) 

5 CI2NTIME 

0 2 7., 1 = 1 ,N 

C0N=P6/ ( XNLJN ( I) ««2«XPN ( I ) ) 

r e a j k = l ,NecE 

P S 1 2 ( I * K ) = PS I ( K * I )**2 
XM2(I ,K) = XMK,I ) ««2 
TAPI C K ) -*- C 2 ( K ) « P S 1 2 ( I * K ) 

TAB2(K)*CCVm*FSI2(I»K) 

CeN5(I,K)»BHe2*FSI2(I»K) 

C £ N6 ( I ,K ) = C£MPS I (K , I ) 

C e N 7 { I , K ) =C£N*XM K » I ) 

C £ N8 ( I ,K )=RH£2*XM2( 1,K ) 

6 = C2MIME 

CALL INTGRA(TA81»N£CE»CX,SUM1( I ) ) 

CALL INTGRA(TAB2,NeCE,CX,SUM2(I)) 

7 CeN TINGE 

De 8-1 I = 1 » N 

C£M (I )= ( F1*XNUN ( I ) «*2«<XMN< I )+P4*SUMl( I ) I » * » 2 

8 • * CEMIME 

IFUFLAG-2) 85,97,85 
85 Cl =6. 198 
C2 = : .84 5 
C3 = :.C:3 
C 4 = 1 . ! 

C 5 =- '.63 
C6 = '1.3886 
DF5 = J.'J36 
GAM* ,‘.7 
GAM2=0.3 
P5 = RHe/(4. ^ * P I ) 

IF(STANT-1.) 197,190,197 
198 WRITE (6 ,6 J52 ) 

6752 F0PMAT(1F1 »4X2JFCEF IN IT I £N 2F 0UTPUT//) 

WRITE (6 ,6335 ) 

6' 35 F2RPAT (4 >4?HC BAR = MEAN DRAG C IS PL ACEH ENT— IN . ) 

WRITE (6,6 ;36) 

6 736 FePNAT (4X49HP BAR = MEAN DRAG 0 ENC ING M0MENT — LB*IN.*E-6) 

WRITE (6 , 6 337 ) 

6:37 FeRRAT ( 4 > 5 2F S IGF A ( C ) = STANCARC DEVIATI0N 0F C ISPLACEPENT— IN. ) 

WRITE (6,6338) 

6738 F2PPAT(4X49FSIGPA(N ) = STANCARC OEVIATI0N 0F BENDING P0PENT , 

UX12H— LE«IN.«E-6) 

WRITE (6,6 39) 

6 ' 3 9 F2RPAT (4>45FVARIANCE(C)= VARIANCE 0F DISPLACEMENT — IN.**2) 

WRITE (6,6041 ) 

6,141 FePMAT(4X39HVARIANCE(M)= VARIANCE 0F BENDING M0MENT , 

11 X 18 H-- (Le«IN.)**2«E-12) 

WRITE (6, 6043 ) 

6043 F2RPAT (4X32HV = WINC SPEED— IN ./SEC ) 

WRITE (6 ,6-049 ) 63 



6 04 9 FOPMAT(4X50HAESC = MX. ABSOLUTE C ISPL ACEMENT ASSOCIATED, 

1 1 X40HW I T F 3 SIGMA LIFT ANC CRAG RESPONSE — IN.) 

WRITE (6,6050) 

6 j 5 0 F OPMAT (4X52HAESM = MAX. ABSOLUTE BENDING MOMENT ASSOCIATED, 

11X47HWITF 3 SIGMA LIFT ANC CRAG RESPONSE— LB • IN .*E-6 ) 

WRITE (6 ,6044) 

6 ''44 FORMAT (4X53HCINT = ABSOLUTE DISPLACEMENT! BENOING MOMENT ) 

11X17H— IN. (LB«IN.«E-6) ) 

WR ITE (6 ,6046 ) 

6H46 FORMAT (4X44HSFANS = PROBABILITY DENSITY 0F ABSOLUTE* 

1 1 X28HD IS PLACEMENT (BENC INC MOMENT ) ) 

WRITE (6,6047) 

6 '47 FORMAT (4X48HFANS * CUMULATIVE DISTRIBUTION OF ABSOLUTE* 

11X28HCISFLACEMENT(EENCINC MOMENT)) 

197 CONTINUE 

ST ANT = ST ANT+1 • 

WRITE (6 ,2010) 

2Clw FORMAT (1H1.25H OUTPUT IN CRAG DIRECTION) 

GO TO 95 
SO C 1 =3 . 954 
C 2 =0. 51 2 
C3 =C • 1 
C4=3 .58 
C5=-3.35 
C6 =0.3886 

DF 5=0. 04 5 .... 

G A Ml = 1 .4 
GAM2-0.5 
P5=RH0/ (8.0-PI) 

WR ITE (6 ,2011) 

2 “1 1 FORMAT (1F1 ,25F OUTPUT IN LIFT DIRECTION) 

95 DO 97 1 = 1, N 

C 2 N3 ( I ) = (XMN(I)/(XMM I)4P4«SUM1( I) ) )* 

* (BETAN (I )+(F5/(XNUN( I )«XMN( I) ) )*SUM2( I ) )**2 
97 CONTINUE 

U = C6#W183 . 

U2=l*U 

C215=1.5«C2 

C253=5.0/(3.0«C2 ) 

UC=U2*C 1 _ 

DO 96 I = 1 , NOCE 

FM(I)=C3«H183(I)«*C4 . . 

BETA(I)=F183(I)««C5 

FF(I)=H(I)/(W(I)*FM(I)) 

96 CONTINUE 

DO 120 J = 1 ,NFREC 
P8 = XNU ( J ) * P3 

P9=F8«P8 .. . .. . 

DO 130 1 = 1, NOCE 

FFM = XNU(J)«FF(I ) 

S < I ) =BET A ( I ) «UC«FFM/(XNU( J )«( 1 .O+C 21 5«FFM**C 2 )**C253 1 

100 CONTINUE .. 

DO 110 11=1, NODE 

DO 110 12=11 , NOCE - - . 

DF=AB$ (XNU( J)«H( 12 ) /W ( 1 2 )-XNU < J ) *M ID/MIID) 

EXFO = EXP (-0.3465«CF/DF5 ) 

IF ( (H( II )+H( 12) )/2 .0-10(3 ) 10 2. 102* 103 

102 GA M = GAM 1 

GO TO 104 

1C3 GAM=GAM2 

104 PHI(I1,I2)=(SCRT(S(I1)«S(I2))«E^P0* 



tees (P2«GAMCF)I«155. .0 
111) C £MI ME 

oe 120 K*1,K 
SUP3(J,K)=0.C 

oe ns i = i,NecE 

SUP3 ( J,K)*SUP3(J,K)+(CCVS( I )4P9«D4< I ) )*PHIH , I )«PSI2(K, I ) 

115 C0KTIME 
SUP4 < J,K >=0.9 

I F ( I FLAG-2 ) 119 » 1 16 » 1 19 

116 D0 117 I = l , M 
11=1+1 

D £ 117 L = 1 1 ♦ MCE 

SIP4 ( J,K )=SUP4< J,K)+(CCV( I) «CCV C L ) /A .0+P9+D2 (I ) «D2f L I l*PHI ( I .1 ) 
S*PSI (I»K)«PSI(L»K) 

117 C2MIME 

G £ T 0 121 

119 D £ 118 1*1, M 

11=1+1 

D £ 118 L= 1 1 » MDE 

SUK( J,K) = SUPA< J,K) + (CCV(I»«CCVai+P9*D2( n«C2U II*PHIil,U. 
«*PSI (I,K)«PSI(L,K) 

118 C0NTIMJE . 

121 C0FTIMJE 

SUP3(J,K)*DX2*ISUP3(J,K)42.0«SUM4<J,K)) 

120 C0MIMJE 

D£ 250 L = i ,N2CE . . . . 

D 2 200 J = 1 »NFPEC 

PHIMJ)=0.0 

PH I Y ( J ) =0 .0 

D £ 140 K = 1 , N 

C0KST= ( SUPS I J »K ) / < C0N 1 < K ) « < < 1 .0- « XNU J ) / 

$XMN f K ) ) *«2)««2 + C2M(K)«4.0MXNU(J) /XNUN CK ))«*2 ) )) 

PHIP(J)=FHIP(J)4C£N8|K,L)«CBNST 

PHIY(J)=FHIY(J)+C£N5IK,L)«C£NST 

140 C0NTIME 

2Cw C 2MI ME 

CALL I NT GRA ( PHI P , NFPE C , C FR EC , V ARM! L )1 

CALI INTGRA(PHIY,MFEC,CFREC,VARC<1 )) 

250 C2NTIME 

IF UFLAG-2)360,380, 360 . 

360 D0 350 1 = 1 * MCE 

SIGCYU )=SCPT(VARC< I)> .... . 

SIGPY(I )=SCRT(VAPP( I!) 

0B(I)*0.0 ... 

XPB ( I ) =Q .0 

00 350 K = l,K 

D 0 275 J=1,K£CE 

275 ARRAY(J)=CCV2(J)*PSI(J,K) 

CALL INTGRAI ARRAY, N£CE»CX,XIM) 

DB U)=CE U) + C£N6(K, II«XINT 

XPB(I) = XPEm + C£K7{K, mXINT 

350 C0MIKUE 

VARC1=VAPC(1) 

VARD2 = VAPC 111 

VARP1 = VAPP ( 1 ) 

VAPP2 = VAPM1) 

DB1 = D8 (1 ) 

DB2 = DB ( 1 ) - 

XPei=XPE(l) 

XHB2=XPE (1) 

00 351 I =2 , N £CE 65 



IF (VaPC { I ) .GT.VARC2 I V A R C2-VAR D( I ) 

IF (VARC (I) .LT.VARCl ) V ARC 1*VARC( I ) 

1 F (VAR MI) .GT.VARP2 5 VAPP2*VAPM(I) 

IF (VAP.M I) . LT.VAPPl ' VAPP1=VARM(I) 

IF (CBU) .GT.CeZ) CB 2 = H E ( II 

IFfCB(I) .LT.Cei) CPI = C E ( I ) 

IF(XPe(I).GT.XME2) X P F 2 =XPE(I) 

IF(XPe(I).LT.XPPl) XPE1 = XP 6 ( I ) 

351 C2NTIMJE 

WRITE (6 ,20 1HCEI I ) t I»1»N2CE ) 

WRITE (6 ,2X»4I (XPei I), I=1,N2CE) 

CALL CUlK3L(-l.Xl,X2,CEl,Ce2,6F B 0 TT, S IDE 1 ,-N0DE , X ,DB ) 

CALL CUlP3l(-l,Xl t X2,XPF 1,XPE2,6F B0TT, S ICE2 »-N0DE » X » XPB ) 

CALL CLIK3L(-1,X1,X2,VAPC1,VAPD2,6H B0 T T , S IDE 3 ,-NBDE , X , VARD ) 

CALL CUIK3U-I.XI v X2»VARP 1,VAPM2,6H , B0TT » $ IDE 5 ,-NBDE , X , VAR P ) 

WRITE (6 ,2000 I ( S IGC V ( i i , 1= 1 ,N0OE ) 

WRITE (6, 2003) ( S IC-P V ( 1 ) , 1= 1 ,N2 DE ) 

G2 T2 39»" 

383 CONTINUE 

02 381 I 3 1 , N 2CE 
SIGCZ(I) 3 SCPT(VARC( I) } 

SIGPZ(I >=SCRT(VAPP( I) ) 

ABSC(I )*XP«SIGDY( I ) 

ABSP ( I )*XR*SIGPY( I ) 

381 CEMIMJE 

VARC1 = VARC ( 1 ) 

VARC2 = VAPC ( 1 ) 

VAPP1 = VA fip ( 1 ) 

VA PP2 = V A P P ( 1 ) 
oe 382 I = 2 * N £CE 

IF ( V AP C ( I ) »GT.VARC2 ) V ARC2=VARC( I) 

IF (VARC (!) .LT.VARCl I VAPC1*V/RC( I) 

IF (VAR P (I ) .GT.VARP2 ) VARP2=VARM(I) 

I F(VARP(I). LT.VAPPl) VARP1=VARM( I) 

382 C0NTIMJE 

WRITE (6,2303 ) ( S I G C Z ( I),I=1,N2CE) 

WRITE (6,2003) (SI GPZ (I ), 1*1, NeCE) 

CALL CUIH3U-1 ,Xl,X2,VARCl,VAPD2,fcH ,B0TT»SIDE4»-N2DE,X,VARD) 

CALL CUIK3L(-1,X1,X2,VARP1,VARM2,6H B0 T T , £ IDE fc ,-N2DE , X , VAR P ) 

390 C2NTIMIE 

WRITE (6,2006) ( V A P C ( I ) , I = 1 , N2 C E ) 

WRITE (6 ,2007) (VARP( I ) , 1= 1 ,N2CE ) 

2006 F2RPAT (50X,11FVARIANCE(C )/( 6E16.8)) 

2'iC7 F£RPAT(5^X,11FVARIANCE(P)/(€E16.0)) 

2000 F2PPAT(5! X f 8FSIGPA(C)/(6E16.E)) 

2C1C1 F 2PP AT ( 5C X , 5FC E AR/ ( 6 E 16 . 5 ) ) 

2C03 F2RPAT (5 3X,8FSIGPA(P)/(6E16«E ) ) 

2004 F2PPAT (5”X,5HP EAR/ (6E 16 .8 ) ) 

WRITE (6 ,1113) (V(I ), I=1,N2CE) 

1113 FBPPAT (/53X,1FV/(6E16«8 ) ) 

IFUB2TF) 393,393,392 
392 CBKTIME 

De 422 I * 1 , N ECE 
C CBPPUTE «*AESP«« 

A*XR*SIGPZ ( I ) 

AB$XPB = AES (XPE ( I ) ) 

IFIA.GT.AESPU)) Ge T2 4 10 
AB$P(I) = ABSXPE+ABSP(I ) 

Ge Te 417 

410 AB2*ABSXPB/(A«»2-A6SP(I)««2) 

IF (AB2*AESP( I) .LT.l .0 ) G2 Te 415^^ 
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ABSMn^esxye+JesFU > 

ce T 0 417 

415 ABSMI )*A«(1.0 + AE$XM : «Ae2)««C.5 

c 

C C0NFUTE ««ABSC«* 

417 A= >P*S IGCZ ( I ) 

ABSCB = AB$ (CB ( I ) ) 

IF(A.GT.AESCU)) G 2 T£ 419 
A8SC ( I ) = ABSCe + AeSC { i ) 

G0 T0 422 

419 AB2=ABSCB/(A«'»2-AESC( li*«2) 

IF <AB2«AESC< II.LT.1.0 ) G0 T0 421 
ABSC(I) = AESCe4AeSC{ I) 

Ge T0 422 

421 ABSC(I) = AM1.0 + AE$CE«AE2)««C.5 

422 C0NTIME 
AB$01*AE$C (1 ) 

AB$C2 S AP$C ( 1 ) 

A6SPl*AeSP (1 ) 

ABSR2 = AB$M1) 

00 423 I=2,NeC£ 

IF (ABSC (T ) .GT.AESC2 ) AESC2=AESC(I) 

IF (A8SCm.LT. AES Cl ) AESC1=AESC(I) 

IF(ABSMI),GT.AESP2) ABSP2*AE$M( I ) 

IF(AB$MI).lT.AESn ) A i = AESM( I ) 

423 C0MIMJE 
WRITE (6 , 6053 ) 

6053 F0PPAT ( 1 FI *4X18F RESULT AM RESF0NSE//) 

WRITE (6, 2021) ( AESC (I ) , I * 1 , N0CE ) 

WRITE (6» 2022) <AB$M( I )» I=1»N0CE) 

2^21 F0PRAT(53X,4HABSC/(6E16.8 ) ) 

2022 F0PMT(5OX »4HABSR/( 6E16.8) ) 

CAU CUIK3L(-1,X1.X2* AES Cl, AESD2,6H B0TT, $IDE7,-N0DE *X,AB$D) 

CALL CUIK3L(-1,X1,X2,AESK1, A ESM2 , 6H. ...... B0TT , S IDE8 ,-N0DE »X »AB SP) 

353 C0KTIMJE 

IFIIFLAC— 3)1,4CK,1 
4C0 I F LAG = 2 
I B 0TH= 1 
G 0 T 0 90 

EVALUATE F AS A FUNCTI0N 0F C ANC M ABS0LUTE 

500 00 800 L*I »2 

IF (L.EC.l) Ge T 0 722 
NI FP-M FP 
G0 T 0 723 

722 NI FF=NI FC 

723 00 800 KS*1»MPF 
IF (L.EC.l ) GET2 72 j 

c C0RPUTE F e R ( F ) 

NIF=IPTMKS) 

SY*SIGKY(MP) . „ 

SZ«SIGPZ(MP) 

SS=SIGRS(KS) . 

SB=SBARMKS) 

ABSCP*ABS(XPE(NIF) ) _ 

0DIM(2)-l«0«AeSP(MP)/KP0 

G0 T 0 430 

C C0PFUTE F2R (C) 

720 NIP*IPTC (KSI . . 
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SZ*SIGCZ IMF) 

SS-SIGCS(KS) 

SB*SBARC <KS) 

AB$CP*AeS(Ce(NIF) ) 

DDINT(1)=1.0«ABSC(NIP)/NFC 
430 CONTINUE 
NCT*f) 

DB 600 1*1, NFC 
NlCT*NCT+l 
SFANS(I)*C.O 
DINT(I)=(I-1)«0CINT(L) 

SS *2*0.0 

SS*4*0.C 

DB 70 0 J = 1 » N FT 

THETA* (J-l) «CTHETA 

Y*CINT(I)«CeS(TFETA) 

Z*CINT (I)«SIN(TPETA ) 

Q*C.5«(< CY-ABS0M/SY)««2-KZ-SE )*»2/( SZ*«2-»SS»*2 ) ) 
XINT*EXF(-C)«CINT(I )/ ( P 2 «SY « ( SZ «*2+$S«« 2 )•*€ . 5 ) 

IF (J.EC.l.BR.J.EC.NFT ) C£ T £ 480 
IFHJ/2) -2.EC.J) GB T £ 460 
SS*2*SS*2+XINT 
G£ T £ 700 

460 SS *4*SS N4 + X I NT 
G £ T £ 700 

480 SPANS (I ) * XINT4SFANSU) 

700 CONTINUE 

SPANS (I )*<OTHET A/3.3 )« (SPANS! I )+4.0«SSM4*2.C*SSM2) 

SPANS (I ) *2«0*SFANS ( I ) 

IF(I.EC.l) G£ T £ 560 

FANSm»<0.5«<SFANS( 1-1 ) + SFANS< 1 1 >«0C IN T ( L )+FANS( 1-1) 

G£ T £ 6C0 
560 FANS(1)*^.0 
600 CONTINUE 

01 NTI*CI NT € 1 ) 

DI NT2*C I NT ( 1 ) 

FANS1*F ANS ( 1 ) 

FANS2*FAN$(1) 

SFANS1*SFANS ( 1 ) 

SF ANS2*SFANS 1 1 ) 

DB 605 1*2, NCT 

IF(CINT1.GT.CINT(I)) CINT1-CINTI I) 

IF (Cl NT2.1T. CINT( I) ) CINT2*CINT( I) 

IFCFANSl.GT.FANSd) ) FANS1*FANS( I ) 

IF(FANS2.LT.FANS(I) ) F ANS2*f ANS( I) 

IF (SPANS1.GT.SFANS( I) ) S F ANS 1 =SF ANS ( I) 

IP(SFANS2.tT.SPANS(I) ) S F ANS2*SFANS( I) 

605 CBNTINUE 

CALL CUIK3L(-1»CINT1»CINT2»FANS1»FANS2»6H ,80TTP(1,L) , 

SSIDE9, -NCT, CINT, FANS) 

CALL CUIK3L(-1,CINT1,CINT2, SPANS 1,SFANS2,6H ,B0TTP( 1,1 > , 

SSI CE10, -NCT, CINT, SPANS) 

MR ITE (6 , 2005 ) ( C INT ( I) , I* 1 , NCT ) 

WRITE (6,2002) (SF ANS ( I ) , I* 1 , NCT ) 

WRITE (6,2318) ( FANS ( I ) , I* 1 , NCT ) 

2005 F0R*AT151X»4HCINT/ (10E13.5)) 

2002 F0R*AT(51X,5PSFANS/( 10E13.5) ) 

2318 F0RRAT (5 IX ,4HFANS/ ( 10E 13 .5 ) ) 

800 CONTINUE 
G0 10 4 
813 CBNTINUE 
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CALL ENCJ2B 
STCF 

END 

FCR,IS IMG, IMG 

SURP01TUE IMGRA(Y,N,CX,ANS) 
DlNFKSm Y(l) 

SUYO.O 

d e n 1 = 1 ,n 

1 ■ > SL> = SIMY(I) 

ANS=DX«SLR 

RETLRf^ 

ENG 

XC7 ,D 
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